﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;

namespace CPAIS.CSFsAnalysis
{
    public class FuzzyDEMATEL
    {
        public Dictionary<int/*expert id*/, double[,]/*evaluation result provided by expert*/> _matricesLow;
        public Dictionary<int/*expert id*/, double[,]/*evaluation result provided by expert*/> _matricesMid;
        public Dictionary<int/*expert id*/, double[,]/*evaluation result provided by expert*/> _matricesUp;
        private double[,] _matrix;
        int _nExpertNum;
        int _FactorNum;

        public FuzzyDEMATEL()
        {
            
        }

        public void analysis(string szDirectoryPath)
        {
            readMatrices(szDirectoryPath);
            _matrix = measureMatrix();
        }

        public void readMatrices(string szDirectoryPath)
        {
            string[] dirs = Directory.GetFiles(@szDirectoryPath, "c*");
            _nExpertNum = dirs.Length;
            for (int nIndex = 0; nIndex < _nExpertNum; nIndex++)
            {
                string szFilePath = dirs[nIndex];
                string[] readText = File.ReadAllLines(szFilePath);
                _FactorNum = readText.Length;
                double[,] itemValuesLow = new double[_FactorNum, _FactorNum];
                double[,] itemValuesMid = new double[_FactorNum, _FactorNum];
                double[,] itemValuesUp = new double[_FactorNum, _FactorNum];
                for (int i = 0; i < readText.Length; i++ )
                {
                    char[] chars = new char[1];
                    chars[0] = ',';
                    string[] szItemsIneachLine = readText[i].Split(chars);

                    for (int j = 0; j < szItemsIneachLine.Length; j++)
                    {
                        string szValue = szItemsIneachLine[j];
                        if(szValue.Equals("NO"))
                        {
                            itemValuesLow[i,j] = 0.0;
                            itemValuesMid[i,j] = 0.0;
                            itemValuesUp[i,j] = 0.25;
                        }
                        else if (szValue.Equals("VL"))
                        {
                            itemValuesLow[i, j] = 0.0;
                            itemValuesMid[i, j] = 0.25;
                            itemValuesUp[i, j] = 0.5;
                        }
                        else if (szValue.Equals("L"))
                        {
                            itemValuesLow[i, j] = 0.25;
                            itemValuesMid[i, j] = 0.5;
                            itemValuesUp[i, j] = 0.75;
                        }
                        else if (szValue.Equals("H"))
                        {
                            itemValuesLow[i, j] = 0.5;
                            itemValuesMid[i, j] = 0.75;
                            itemValuesUp[i, j] = 1.0;
                        }
                         else if (szValue.Equals("VH"))
                        {
                            itemValuesLow[i, j] = 0.75;
                            itemValuesMid[i, j] = 1.0;
                            itemValuesUp[i, j] = 1.0;
                        }
                    }
                      
                }
                _matricesLow.Add(nIndex, itemValuesLow);
                _matricesMid.Add(nIndex, itemValuesMid);
                _matricesUp.Add(nIndex, itemValuesUp);
            }

        }

        public double[,] measureMatrix()
        {
            double[,] minLow = extractMinLow();
            double[,] maxUp = extractMaxUp();
            Dictionary<int, double[,]> matrices = new Dictionary<int, double[,]>();
            for(int nIndex = 0; nIndex < _nExpertNum; nIndex++)
            {
                double[,] crispValues = new double[_FactorNum, _FactorNum];

                for (int i = 0; i < _FactorNum; i++)
                {
                    for (int j = 0; j < _FactorNum; j++)
                    {
                        double deltaMinMax = maxUp[i, j] - minLow[i, j];
                        double lambda = 0.0;
                        double mu = 0.0;
                        double rho = 0.0;

                        if (deltaMinMax > 0)
                        {
                            lambda = (_matricesLow[nIndex][i, j] - minLow[i, j]) / deltaMinMax;
                            mu = (_matricesMid[nIndex][i, j] - minLow[i, j]) / deltaMinMax;
                            rho = (_matricesUp[nIndex][i, j] - minLow[i, j]) / deltaMinMax;
                        }

                        double omega = mu / (1 + mu - lambda);
                        double varphi = rho / (1 + rho - mu);
                        double x = (omega * (1 - omega) + varphi * varphi) / (1 - omega + varphi);
                        double z = minLow[i, j] + x * deltaMinMax;
                        crispValues[i, j] = z;
                    }
                }
                matrices.Add(nIndex, crispValues);
            }

            double[,] ret = new double[_FactorNum, _FactorNum];

            for (int i = 0; i < _FactorNum; i++)
            {
                for (int j = 0; j < _FactorNum; j++)
                {
                    double dTotal = 0.0;
                    foreach (double[,] matrix in matrices.Values)
                    {
                        dTotal += matrix[i, j];
                    }
                    dTotal = dTotal / _nExpertNum;
                    ret[i, j] = dTotal;
                }
            }

            //step 2. normalize
            double dMax = double.MinValue;
            for (int i = 0; i < _FactorNum; i++)
            {
                double dRowTtotal = 0.0;
                for (int j = 0; j < _FactorNum; j++)
                    dRowTtotal += ret[i, j];
                if (dRowTtotal > dMax)
                    dMax = dRowTtotal;
            }
            for (int i = 0; i < _FactorNum; i++)
            {
                for (int j = 0; j < _FactorNum; j++)
                    ret[i, j] = ret[i, j] / dMax;
            }

            //step 3. Constructing the total relation matrix (T) 
            double[,] identifyMatrix = new double[_FactorNum, _FactorNum];
            for (int i = 0; i < _FactorNum; i++)
            {
                for (int j = 0; j < _FactorNum; j++)
                {
                    if (i == j)
                        identifyMatrix[i, j] = 1;
                    else
                        identifyMatrix[i, j] = 0;
                }
            }
            double[,] subRet = this.sub(identifyMatrix, ret);
            double[,] inverseRet = this.inverse(subRet);
            double[,] multipleRet = this.multiple(ret, inverseRet);
            return multipleRet;
        }

        public double calR(int i)
        {
            if (_matrix == null)
                return double.NaN;
            else
            {
                double dret = 0.0;
                for (int j = 0; j < _FactorNum; j++)
                    dret += _matrix[i, j];
                return dret;
            }
        }
        public double calD(int j)
        {
            if (_matrix == null)
                return double.NaN;
            else
            {
                double dret = 0.0;
                for (int i = 0; j < _FactorNum; j++)
                    dret += _matrix[i, j];
                return dret;
            }
        }
        public double calP(int i)
        {
            return calR(i) + calD(i);
        }
        public double[] calP()
        {
            double[] ret = new double[_FactorNum];
            for (int i = 0; i < _FactorNum; i++)
                ret[i] = calP(i);
            return ret;
        }
        public double calE(int i)
        {
            return calR(i) - calD(i);
        }
        public double[] calE()
        {
            double[] ret = new double[_FactorNum];
            for (int i = 0; i < _FactorNum; i++)
                ret[i] = calE(i);
            return ret;
        }

        public double calThresholdValue()
        {
            double total = 0; 
            for (int i = 0; i < _FactorNum; i++)
                for(int j = 0; j < _FactorNum; j++)
                total += _matrix[i,j];
            double average = total / _FactorNum * _FactorNum; //平均数             
            double fto = 0;
            for (int i = 0; i < _FactorNum; i++)
                for(int j = 0; j < _FactorNum; j++)
                    fto += Math.Pow((average - _matrix[i, j]), 2);
            double std = fto / _FactorNum * _FactorNum; //方差
            return average+std;
        }
        private double[,] extractMinLow()
        {
            double[,] dRet = new double[_FactorNum, _FactorNum];
            for (int i = 0; i < _FactorNum; i++)
            {
                for (int j = 0; j < _FactorNum; j++)
                {
                    double dMin = double.MaxValue;
                    foreach (double[,] valueItems in _matricesLow.Values)
                    {
                        if (valueItems[i, j] < dMin)
                            dMin = valueItems[i, j];
                         
                    }
                    dRet[i, j] = dMin;
                }
            }
            return dRet;
        }
        private double[,] extractMaxUp()
        {
            double[,] dRet = new double[_FactorNum, _FactorNum];
            for (int i = 0; i < _FactorNum; i++)
            {
                for (int j = 0; j < _FactorNum; j++)
                {
                    double dMax = double.MinValue;
                    foreach (double[,] valueItems in _matricesLow.Values)
                    {
                        if (valueItems[i, j] > dMax)
                            dMax = valueItems[i, j];

                    }
                    dRet[i, j] = dMax;
                }
            }
            return dRet;
        }

#region"Matrix operators"
         private double[,] sub(double[,] a, double[,] b)
         {

             double[,] c = new double[_FactorNum, _FactorNum];
             for (int i = 0; i < _FactorNum; i++)
             {
                 for (int j = 0; j < _FactorNum; j++)
                 {
                     c[i, j] = a[i, j] - b[i, j];
                 }
             }
             return c;
         }

          private double[,] multiple(double[,]  a, double[,]  b)
          {
              double[,] c = new double[_FactorNum, _FactorNum];
              for (int i = 0; i < _FactorNum; i++)
              {
                  for (int j = 0; j < _FactorNum; j++)
                  {
                      c[i, j] = 0;
                      for (int k = 0; k < _FactorNum; k++)
                      {
                          c[i, j] += a[i, k] * b[k, j];
                      }
                  }
              }
              return c;
          }
          private double[,] inverse(double[,] m)
          {
             //clone
             double[,] a = new double[_FactorNum, _FactorNum];
             for (int i = 0; i < _FactorNum; i++)
             {
                 for (int j = 0; j < _FactorNum; j++)
                 {
                     a[i, j] = m[i, j];
                 }
             }
             double[,] c = new double[_FactorNum, _FactorNum];
             for (int i = 0; i < _FactorNum; i++)
             {
                 for (int j = 0; j < _FactorNum; j++)
                 {
                     if (i == j)
                         c[i, j] = 1; 
                     else 
                         c[i, j] = 0; 
                 }
             }
             
             //i表示第几行，j表示第几列
             for (int j = 0; j < _FactorNum; j++)
             {
                 bool flag=false;
                 for (int i = j; i < _FactorNum; i++)
                 {
                     if(a[i,j]!=0)
                     {
                         flag = true;
                         double temp;
                         //交换i,j,两行
                         if (i != j)
                         {
                             for (int k = 0; k < _FactorNum; k++)
                             {
                                 temp = a[j, k];
                                 a[j, k] = a[i, k];
                                 a[i, k] = temp;
 
                                 temp = c[j, k];
                                 c[j, k] = c[i, k];
                                 c[i, k] = temp;
                             }
                         }
                         //第j行标准化
                         double d = a[j, j];
                         for (int k = 0; k < _FactorNum; k++)
                         {
                             a[j,k]=a[j,k]/d;
                             c[j, k]= c[j,k]/d;
                         }
                        //消去其他行的第j列
                         d = a[j, j];
                         for (int k = 0; k < _FactorNum; k++)
                         {
                             if(k!=j)
                             {
                                 double t = a[k, j];
                                 for (int n = 0; n < _FactorNum; n++)
                                 {
                                     a[k,n]-=(t/d)*a[j,n];
                                     c[k,n]-=(t/d)*c[j,n];
                                 }
                             }
                         }
                     }
                 }
                 if (!flag) 
                     return null;
             }
             return c;
         }

#endregion
    }
}
